Topological constraints on spiral wave dynamics in spherical geometries with 

inhomogeneous excitability 



O 

o 

(N 

-t— > 
o 

O 

oo 



Jorn Davidsen/' -^'0 Leon Glass, '^'Q and Raymond KapraP'^'0 

' Max- Planck- InsUtut fiir Physik Komplexer Systeme, 
Nothnitzer Strasse 38, 01187 Dresden, Germany 
^Chemical Physics Theory Group, Department of Ghemistry, 
University of Toronto, Toronto, ON M5S 3H6, Ganada 
"'Department of Physiology, McGill University, Montreal, Quebec H3G 1 Y6, Ganada 

(Dated: 8th February 2008) 

We analyze the way topological constraints and inhomogeneity in the excitability influence the 
dynamics of spiral waves on spheres and punctured spheres of excitable media. We generalize the 
deflnition of an index such that it characterizes not only each spiral but also each hole in punctured, 
oriented, compact, two-dimensional differentiable manifolds and show that the sum of the indices is 
conserved and zero. We also show that heterogeneity and geometry are responsible for the formation 
of various spiral wave attractors, in particular, pairs of spirals in which one spiral acts as a source 
and a second as a sink — the latter similar to an antispiral. The results provide a basis for the 
analysis of the propagation of waves in heterogeneous excitable media in physical and biological 
systems. 

PACS numbers: 87.19-j, 05.40.-a, 89.75.-k 
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I. INTRODUCTION 

Geometry and inhomogeneity influence pattern forma- 
tion in chemical and biological systems. P, One ex- 
ample where these two factors play a crucial role is in 
the experimental observations of distinctive spiral wave 
dynamics on the surfaces of spherical beads, which are 
excitable inhomogeneous chemical media A biolog- 

ical example is the origin of abnormal cardiac rhythms in 
the human heart which depend on the anatomical sub- 
strate. The heart possesses a complex nonplanar geom- 
etry with multiple chambers, with holes corresponding 
to valves and blood vessels. Some serious arrhythmias 
are associated with circulating spiral waves similar to 
those observed in chemical media i5]. Since an abnor- 
mal anatomical substrate is a common finding in patients 
with some types of cardiac arrhythmias, and interven- 
tions that_ modify the anatomy are an accepted form of 
therapy 6] , theoretical analyses of the interplay between 
geometry of the substrate and dynamics may help in the 
therapy of cardiac arrhythmias. 

In this paper we study spiral wave dynamics on (punc- 
tured) spheres with spatially inhomogeneous excitability. 
We show for punctured spheres that the sum of indices 
which characterize each spiral has to be zero. Moreover, 
we demonstrate that topological constraints imposed by 
the spherical geometry and inhomogeneity in excitability 
lead to the formation of pairs of spirals, with distinc- 
tive transient dynamics or as stable attractors, in which 
one spiral acts as a source and a second as a sink lead- 



ing to a source-sink pair under a broad range of condi- 
tions. Our results explain the experimental observations 
of spirals on spherical beads 3, 4] . While we do not con- 
sider detailed models of cardiac wave propagation, our 
results may apply to some generic aspects of atrial ar- 
rhythmias because the thin walls of the atria can be de- 
scribed as two-dimensional (2d) inhomogeneous excitable 
media with specific geometrical features. ^ 



II. AN INDEX THEOREM FOR PHASE 
SINGULARITIES 

The mathematical description of spiral waves is based 
on the notion of phase which in turn allows one to char- 
acterize spiral waves by an index. From this description, 
a number of topological results placing restrictions on 
spiral wave dynamics can be derived. IM IS 121, llCj, |ll| 

With the exception of a finite number of singular 
points, with each point in an orientable and compact 
two-dimensional differentiable manifold M we identify a 
unique phase lying on the unit circle, $ G S*^. The re- 
sulting phase map or phase field is assumed to be contin- 
uously differentiable, except at the singular points. The 
manifold can be triangulated |12| | (subdividing it into a 
set of polygons), where none of the edges or vertices of 
the polygons pass through a singularity. The index, / 
(sometimes also called the topological charge or winding 
number) of a curve C bounding a polygon is found by 
computing the line integral 
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where the polygon is always traversed in a clockwise ori- 
entation. By continuity of V<i>, / must be an integer. 
The index of a singular point is uniquely defined as the 
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index of any curve C provided that C encircles the point 
but no other singular points. The index of a curve that 
does not enclose any singular points is obviously zero. 

If the manifold M has no boundaries, each edge of the 
triangulation is an edge of two polygons. Since the line 
integral adds up the change in phase along the various 
edges of the polygon, the sum of the indices of the sin- 
gular points for a phase field in AI is 



Sm/ = 0, 



(2) 



where the sum is over all the singular points. This fol- 
lows since the contribution of the change in phase of each 
edge to the total integral is counted twice, but since the 
edge is traversed in opposite directions each time, the net 
contribution of each edge is zero. This index theorem is 
applicable to tori, and other surfaces of genus different 
from 0. However, unlike the more familiar Poincare index 
theorem (see Ref. ^3 P- '^^ foi' vector fields) the sum of 
the indices of the singular points does not depend on the 
genus of the surface. 

This index theorem for manifolds without boundaries 
can be extended to manifolds with boundaries. In the fol- 
lowing, we will consider the case of structures that arise 
from puncturing orientable and compact two-dimensional 
differentiable manifolds. The index of a hole can be 
uniquely defined as the index of a curve C provided that 
C encircles the hole but no other holes or singular points 
and C is positively oriented with respect to the domain 
which contains the hole and is bounded by C . Apply- 
ing this definition and taking the summation in Eq. (0) 
over the singularities and the hole, or the holes if there is 
more than one hole, the index theorem can be proven by 
the same line of arguments as for the case of manifolds 
without boundaries. 

This extension is important in the heart, for example, 
where the atrium is punctured by valves and veins. In 
such cases one is led to consider manifolds with holes; for 
example, a sphere with a hole. A sphere with a hole is 
topologically equivalent to a disk, and, indeed, the results 
for disks and for spheres with holes are consistent: For 
the disk D^, bounded by a curve C, 
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so that the sum of the indices of the singular points in 
the disk is equal to the index of the curve C bounding 
the disk. If there is a single singular point on the disk, 
with an index of -1-1, the index of the curve bounding the 
disk will also be -1-1. Imagine now the boundary of the 
disk to be brought together (like a draw-string bag) so 
that the boundary of the disk now defines a hole in the 
sphere. In this geometry the curve C will be traversed in 
an opposite orientation (the hole is now inside C) from 
the direction it was traversed when it was the boundary 
of the disk. Now if there is a singular point with an index 
of -1-1 on the sphere, the index of the hole is -1, so that 
the sum of the indices is again zero. 



Since it is necessary to conserve the sum of the indices, 
singularities of index ±1 usually arise and are destroyed 
in pairs of opposite sign . An exception occurs when 
singularities are destroyed by collision into a boundary, so 
that the index of the singular point and the index of the 
boundary both change simultaneously. Another excep- 
tion occurs if there are singularities with index different 
from one. In such cases interactions between different 
singularities can lead to destruction or creation of odd 
numbers of singularities |16|. 



III. THE FITZHUGH-NAGUMO EQUATION 



The FitzHugh-Nagumo (FHN) equation, [T] 



du 



(4) 



where u{r,t) and v{r,t) are two scalar fields, is the 
ratio of the time scales associated with the two fields, 
and I?„ and are the constant diffusion coefficients, 
is a prototypical model for excitable media. We choose 
Du — 2 and Dy = 0. The real parameters a and 
(3 characterize the local dynamics and, hence, the lo- 
cal excitability. Whenever < a < 1, ae^ < 1 and 
\P\ > Ph = {l- ae2)i/2(i/3(2Q, + a^e^) - 1), the FHN 
system is excitable. At I3h a Hopf bifurcation occurs 
such that for < the system exhibits oscilla- 
tions. In the following we take a = 0.2, e = 0.2 and 
(3> (3h = 0.863.... 

We consider a spherical shell whose outer and inner 
radii are Re and Ri , respectively, and focus on thin spher- 
ical shells where Ri = 40, Re = 42. The radii are large 
enough to avoid a curvature-dependent loss of excitability 
[l^ , and the shell is sufficiently thin that the dynamics is 
effectively 2d corresponding to the dynamics on a sphere. 
[T^ The initial condition is a domain of "excited" state, 
adjacent to a domain of the "refractory" state. Both do- 
mains have the forms of slices of the same size oriented 
from the north to the south pole and yield a pair of 
counter-rotating spirals. 

In order to apply the topological results to the FHN 
equation, it is necessary to first define the phase. We de- 
fine a phase, <I>(r,t), based on the equation tan<I>(r,t) — 
v{r,t)/u{r,t) if v{r,t) ^ and u{r,t) ^ 0. Thus, sin- 
gular points at given t are points r in the medium for 
which v{r, t) — and u{r, t) = 0. For each t, we obtain 
a continuously differentiable phase map AI* — $(:,t)|Dt 
that associates to each point in a well-defined domain 
2?* a phase lying on the unit circle, $65^. For our 
FHN medium, the domain is the surface of a (punctured) 
sphere reduced by a finite number of points where the 
phase is singular at fixed t. 

Rotating spiral waves in the FHN equations are ob- 
viously associated with a singular point which is called 



3 




N 




Figure 1; (Color online) Left: Spiral waves of excitation 
(light fronts) on the sphere for constant /3 — (3ex emanat- 
ing from spiral cores close to the poles on an equator pro- 
jection. The white arrows show the direction of propagation. 
One annihilation front along the equator can be identified. 
Right: Sketch of the constant gradient in the inhomogeneous 
case. The dashed lines are the equi-/9 lines and we choose 
Pmax — 1.0 and Pmin = Pex ■ The angle 4> describes the orien- 
tation with respect to the axis from pole to pole. The results 
described in the text do not depend qualitatively on the choice 
of (p or fUmin and /3max as long as they yield stable spirals. 

the spiral core. In what follows, we assume that there 
are only single-armed spirals so that a clockwise rotating 
source has an index of -1-1, and a counterclockwise rotat- 
ing source has index -1. A clockwise rotating sink has 
an index of -1, and a counterclockwise rotating sink has 
index -1-1. From Eq. (|2J), it is impossible to have a single 
rotating spiral wave on a sphere; in addition, there must 
be at least another singular point or a hole with nonzero 
index. 

For excitable media, a non-zero index of a hole implies 
that wave fronts travel permanently around the hole such 
that the numbers of fronts travelling clockwise and coun- 
terclockwise are different. This includes the particularly 
simple case of a single wave front travelling around the 
hole which can be considered as a spiral wave associated 
with the hole. 



IV. SPIRAL WAVE DYNAMICS IN SPHERICAL 
GEOMETRIES 

A. Dynamics in excitability gradients 

First consider homogeneous FHN media with a con- 
stant P = Pex = 0.9. The wave fronts emanating from 
the spiral cores with opposite index "annihilate" along 
the equator such that each spiral determines the dynam- 
ics on half of the sphere (see Fig.^ leftjpanel) — similar 
to what has been observed in Ref. 21] for a different 
excitable system. We have shown that this behavior is 
robust with respect to disorder in excitability with small 
amplitude and correlation length. If random, uncorre- 
lated spatial variations in [3 exist on length scales much 
smaller than the diameter of the spiral core meander , 
the dynamics is able to average over such small-scale in- 
homogeneities. The robustness explains why such states 
have been experimentally observed in some chemical re- 
actions on spherical beads which are intrinsically inho- 
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Figure 2: (Color online) The local excitability f3{r{t)) at 
the spiral cores versus time. Gradient-induced motion of the 
two spiral cores leads to a change in the local excitability at 
the cores with time. The spiral period in the final state is 
To — 13.2 ± 0.1. Four difi^erent regimes can be identified (see 
text). Inset: The final bound pair of counter-rotating spirals 
in regime III for (f) = 51.0° is shown on an equator projection 
such that the point of lowest excitability lies on the central 
longitude. The spiral closer to the equator has index —1 while 
the other one has index -|-1. 



mogeneous [JJ. 

Applying a constant gradient in /3 as sketched in the 
right panel of Fig. leads to a different scenario. The 
time evolution of the spiral pair may be partitioned into 
four distinct regimes as shown in Figs. [3 and O Because 
of the gradient, the frequencies of the two spirals differ 
since a higher value of (3 corresponds to lower excitabil- 
ity, which generally implies a lower spiral frequency |23j| . 
During a short transient T, the spiral with the higher 
frequency assumes control of the dynamics on the 
sphere. The location on the sphere, where wave fronts 
emanating from the two spiral cores annihilate, moves 
toward the core of the low-frequency spiral. Finally, the 
low-frequency spiral core is pushed farther from the high- 
frequency spiral core |24L l25l | (see Fig.|3J). After this short 
transient, the wave fronts travel from pole to pole leading 
to the creation of a source-sink pair. This (intermediate) 
state is shown in Fig. 0] and corresponds to regime I in 
Figs. 121 and 121 Viewed from the low excitability end of 
the sphere, the waves wind into a small region about the 
core, reminiscent of the structure of antispirals, i.e., in- 
ward moving spirals seen in oscillatory media |26| . How- 
ever, the origin of this inward spiral motion in oscillatory 
media differs and is distinct from that observed here. In 
oscillatory media, either spirals or antispirals are stable 
depending on system parameters and the wave length 
diverges on the border in the parameter space between 
these two regimes. Thus, antispirals exist independently 
of spirals. This is not the case here because the genera- 





Figure 3: (Color online) The distance between the two spiral 
cores d{t) versus t. Gradient-induced motion of the two spiral 
cores leads to a change in the distance d between the cores 
with time. The lower and upper curves correspond to the 
distance in and S^, respectively. The dashed lines are the 
respective upper bounds given by the size of the sphere. 



tion of an inward moving spiral relies on the presence of a 
spiral source and spherical geometry. For example, con- 
sider the FHN system with a disc geometry and a radial 
gradient in excitability such that the maximum value of /3 
is located in the center of the disc. In this case, a source- 
sink pair cannot occur because the high-frequency spiral, 
acting as a source, would push the other (low-frequency) 
spiral out of the system, excluding the presence of any 
strong random inhomogeneities in excitability which may 
pin the low-frequency spiral and prevent its motion (see, 
e.g., 01). The lower panel of Fig. ^ shows that rem- 
nants of the low frequency spiral persist in a small area 
around the core. It has been speculated that antispiral 
waves might occur in cardiac tissue 28] . While excitable 
media cannot support a regime of exclusive antispirals 
due to their excitable character, our results show that 
source-sink pairs with similar characteristics could form 
in the heart where the underlying dynamics is excitable, 
the medium is inhomogeneous and the topology is similar 
to a (punctured) sphere. 

Spiral dynamics of the type described above has been 
observed by Maselko and Showalter 0] in experiments on 
the excitable Belouzov-Zhabotinsky reaction on spherical 
beads. They attributed the generation of spiral source- 
sink pairs to inhomogeneities in the medium related to 
differing chemical environments. This is consistent with 
our findings for systems with gradients and is further 
confirmed by the work reported in Ref. 21] . While the 
generation of source-sink pairs due to a gradient in the 
FHN medium investigated here is only an intermediate 
state (but one that persists for approximately 240 spiral 
periods in our simulations), random spatial variations of 
the excitability with a correlation length comparable to 



Figure 4: (Color online) Waves of excitation on the sphere in 
regime I of the gradient-induced dynamics shown in Figs. |21 
and El A source-sink pair has formed. For random spatial 
variations of excitability with a correlation length comparable 
to the diameter of the spiral core meander, the final state is 
very similar to the one shown here Upper panel from 

left to right: View centered at the north pole, south pole and 
the equator. The source at the south pole has index —1 and 
the sink at the north pole index -1-1. The black circles show 
possible choices of C. The white arrow shows the direction 
of wave propagation. Lower panel: Dynamics at the north 
pole. Time increases from left to right with At = 2.5 between 
snapshots. 



the diameter of the spiral core meander or larger can lead 
to a final state consisting of such a source-sink pair . 
This is due to the fact that the source can be trapped in 
a region of depressed local excitability. 

Not only does the gradient in the FHN medium change 
the local excitability but it also induces a drift of the spi- 
ral cores j29j. For our model, the drift is rather slow 
compared to the transition to the source-sink pair which 
takes place during the transient regime T. This can be 
seen in Fig. El (The fluctuations in /3(r(i)) are due to 
spiral meandering.) In regime I, the dominating spiral 
drifts toward lower excitability and its wave fronts con- 
tinuously push the other core in the opposite direction, 
thus, keeping the distance d between the cores approxi- 
mately constant (see Fig.EJ. The fluctuations in d{t) are 
again due to meandering of the spiral cores. Because of 
this drift, the local excitabilities at the two spiral cores 
approach each other until they become equal. 

At this point, regime II in Figs. (21 and is entered. The 
dynamics change drastically: the slaved spiral reverses its 
drift direction and regains control over its own dynamics. 
Both spirals drift toward lower excitability. Due to the 
geometric constraints imposed by the spherical geometry, 
the spirals approach each other until they form a bound 
pair (at t « 7500). 

They finally reach a stable state (at t « 9500) corre- 
sponding to regime III in Figs. [3 and Neither the av- 
erage distance between the spirals nor the average local 



excitability changes further. Yet, on top of the persisting 
unsynchronized meandering of the two spiral cores, the 
bound pair slowly moves along a (closed) equi-/3 curve 
on the surface of the sphere. The direction of the motion 
depends on the initial condition, i.e., whether the spiral 
with positive index was initially closer to the point of 
lowest excitability or to the point of highest excitability 
than the spiral with negative index. The wave dynam- 
ics generated by this bound pair is shown in the inset of 

Fig.ia 

While kinematic theory applies only to spirals with 
large cores, it is instructive to note that this theory pre- 
dicts that the direction of the drift due to gradients de- 
pends on the model system and its parameters Al- 
though our simulations have shown the same direction of 
the drift for a range of parameters in the meander region 
of the FHN phase diagram, it is conceivable that, under 
different circumstances favoring a drift toward higher ex- 
citability, one spiral could act as a permanent source and 
a source-sink pair could be the final attractor. Such a 
scenario would also be consistent with the experimental 
findings in Ref. 0- 
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Figure 5: (Color online) Waves of excitation on the punctured 
sphere propagating towards the hole. A source-sink spiral 
wave pair has formed in the final state. Upper panel from 
left to right: View centered at the north pole, south pole and 
the equator. Lower panel: Dynamics at the south pole. Time 
increases from left to right with At = 2.5 between snapshots. 



B. Punctured spheres 

Next, we consider a homogeneously excitable sphere 
with a single hole. Two scenarios can be observed de- 
pending on the location of the hole with respect to the 
spiral pair. If a spiral wave is not permanently attached 
to the hole, the dynamics is very similar to the case with- 
out any hole. If one of the spirals is permanently attached 
to the hole, the frequency of this spiral is lowered. The 
size of the hole determines the frequency of the spiral be- 
cause the wave front has to travel around the hole. The 
transient dynamics is similar to that in regime T for the 
case with a gradient; however, no drift of the spiral cores 
is induced and the final state is a spiral source-sink pair 
as shown in Fig. [S] Not only is the net index conserved 
during the transition to a spiral source-sink pair but so is 
the index of the individual spirals: during the transition 
an outgoing counterclockwise (clockwise) oriented spiral 
is converted into an ingoing clockwise (counterclockwise) 
oriented spiral. Thus, the formation of spiral source-sink 
pairs conforms with the topological constraints. 

If a gradient as well as a hole is present, the spiral drift 
discussed earlier also determines the final state, which 
depends on the hole's location in the gradient field. For 
instance, if the point of lowest excitability in the medium 
is on the hole boundary, simulations show that the spiral 
which is not attached to the hole will stabilize close to 
this point. In this final state odd numbers of wavefronts 
are attached to the hole, again conserving the topological 
charge of the hole. 



V. CONCLUDING REMARKS 

Inhomogeneities due to spatially-varying excitability 
on (punctured) spherical shells lead to complex spiral 
wave dynamics and the formation of source-sink spiral 
pairs in excitable media. The results presented here are 
immediately applicable to excitable media in more com- 
plicated geometries such as tori or multi-holed tori and to 
situations in which multi-armed spirals are found. This 
includes mathematical modelling of cardiac tissue. The 
approach taken in this paper stresses constraints and as- 
pects that apply to, and must be observed in, all realistic 
models of the heart satisfying certain criteria of continu- 
ity. There are also implications for the treatment of car- 
diac arrhythmias. In cardiology it is sometimes possible 
to develop maps showing the timing of the excitation over 
limited regions of heart ^J. In this case, a sink might be 
confused for a source (of the arrhythmia), and this might 
have implications for the diagnosis of the mechanism and 
the choice of therapy. The current work shows how par- 
tial knowledge about what is happening in some regions 
that could be observed, might be helpful in establishing 
properties of dynamics that could not be observed. While 
the types of sinks we have described here have only been 
observed in chemical media ^3, 4] so far, we certainly ex- 
pect their existence in the cardiological domain. 

We thank F. Chavez, M. Bar, S. G. Whittington and 
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providing numerical tools. This work was supported in 
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